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Background Of The Invention 

[001] 1. Field of the Invention, 

[002] This invention relates generally to satellite positioning system receiver and more 
practically to altitude aided satellite position receiver. 

[003] 2. Related Art 

[004] Satellite positioning systems (SATPS) have been developed that enable SATPS 
receivers to determine a location upon receiving a SATPS spread spectrum signal from a 
plurality of satellites. An example of a SATPS is the Global Positioning Satellite (GPS) system 
maintained by the United States Government. A SATPS receiver must acquire four satellites in 
order to determine the location (latitude, longitude and altitude) of the SATPS receiver. 
However, in urban and rural indoor or outdoor environments, fewer than the required satellites 
may be acquired. Further, strong signal shading and signal fading may adversely affect a 
SATPS receiver ability to acquire the satellites. If four satellites are acquired, then a three- 
dimensional location may be determined. If three satellites are acquired, then a two-dimensional 
location (latitude and longitude) may be determined if altitude is known. Accurate two- 
dimensional location with three satellites is further impacted if the altitude is not adequately 
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known, for example, when the SATPS receiver is used in a region of uneven elevation, such as 

in the mountains. One approach to assist with location determination if three satellites are 

acquired employs an algorithm that uses the last known altitude calculated from four acquired 

satellites. However, the more the altitude of the SATPS receiver has changed from the last 

known altitude, the greater the error that results in the location of the SATPS receiver location 

determination. 

[005] Therefore, there is a need for methods and systems for determining the location of a 
SATPS receiver when only three satellites are acquirable, that overcomes the many 
disadvantages set forth above and others previously experienced. 



Summary 

[006] Systems consistent with the present invention provide an approach to determine 
location in a SATPS when only three SATPS satellites are acquired with the use of digital 
altitude data. The three SATPS measurements are obtained from SATPS satellites and digital 
altitude data is obtained from a terrain database. Three pseudorange equations and an altitude 
equation are used to determine the location of the SATPS receiver. 

[007] Other systems, methods, features, and advantages of the invention will be or will 
become apparent to one with skill in the art upon examination of the following figures and 
detailed description. It is intended that all such additional systems, methods, features and 
advantages be included within this description, be within the scope of the invention, and be 
protected by the accompanying claims. 
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Brief Description of the Figures 

[008] The components in the figures are not necessarily to scale, emphasis instead being 
placed upon illustrating the principles of the invention. In the figures, like reference numerals 
designate corresponding parts throughout the different views. 

[009] FIG. 1 illustrates a satellite positioning system (SATPS) with four satellites and a 
SATPS receiver. 

[010] FIG. 2 is a block diagram of a SATPS receiver of FIG. 1 in communication with a 
server having digital terrain elevation data. 

[011] FIG. 3 is a diagram of bilinear interpolation employed with digital terrain elevation 
data of FIG. 2. 

[012] FIG. 4 is an illustration of the four orthometric heights identified in FIG. 3. 

[013] FIG. 5 is a flow chart of the steps for determining location with a SATPS receiver 

with digital terrain elevation data. 

Detailed Description 
[014] Turning first to Figure 1, that figure illustrates a satellite position system (SATPS) 
100 with a SATPS receiver 102 located on the Earth 104 in communication with four satellites 
106, 108, 1 10, and 1 12. The SATPS receiver 102 via antenna 1 14 receives the SATPS spread 
spectrum signal 1 16, 1 18, 120, and 122 from the four satellites 106, 108, 110, and 1 12. The four 
SATPS spread spectrum signals 116, 118, 120, and 122 are transmitted at the same frequency. 
It is desirable for the SATPS receiver 102 to be receiving SATPS spread spectrum signals from 
as many satellites as possible, but to resolve the location, such as the x, y and z or latitude, 
longitude and altitude coordinates, at least three SATPS satellites in the SATPS system 100 
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must be in communication with the SATPS receiver 102 along with knowledge of the altitude of 

the receiver. Otherwise at least four SATPS satellites are required. 

[015] As an alternative to a SATPS spread spectrum signal from one of the satellites 

106, 108, 110, or 112, an Earth based pseudollite 126 may transmit a SATPS spread spectrum 
signal 128 to the SATPS receiver 102. The SATPS receiver 102 still requires three spread 
spectrum signals and knowledge of altitude in order to determine the location of the SATPS 
receiver 102. In another implementation, a combined wireless device and SATPS receiver 102, 
such as in a cellular telephone, may receive altitude information from the pseudollite 126. In 
other implementations, the altitude information may be transmitted by a wireless data network 
such as a PCS system, Bluetooth system, private wireless data network, and within a television 
carrier wave to name but a few examples. 

[016] In FIG. 2, a block diagram of a SATPS receiver 102 of FIG. 1 in communication 
with a server 250 having digital terrain elevation data is shown. The SATPS receiver 102 
acquires three SATPS spread spectrum signals, such as 1 16, 1 18 and 120 at the radio transceiver 
202 via antenna 1 14. The SATPS spread spectrum signals are filtered by filter 206 under the 
control of navigation processor 208. The navigation processor 208 is coupled to the filter 206, 
clock 210, memory 212, and interface 216. The clock 210 may provide a plurality of clock 
signals from an oscillator to aid in the acquisition and processing of the SATPS spread spectrum 
signals and additionally provides timing for the navigation processor to transfer data to and from 
the filter 206, memory 212 and interface 216. The interface 216 may be coupled to a display 
214 or the location output data may be outputted by interface 216. 

[017] The terrain database is not communicated from the server 250 to the SATPS receiver 
102 in the current implantation. Rather the code phases are sent from the SATPS receiver 102 
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to the server 250. These code phases are converted to pseudo-range measurements in the server 

250. The computation of position determination using the terrain database is not done in 208. 

Rather the computation is done in a processor of 250, such as controller 256. 

[018] If the SATPS receiver 102 is able to only receive or acquire three SATPS spread 

spectrum signal, such as 116, 118 and 120, the navigation processor 208 provides code phase 

information associated with the SATPS satellite signals to the server 250. The server 250 is 

shown in radio communication with the SATPS receiver 102 via server transceiver 252. The 

server 250 has a controller 256 connected to the server transceiver 252, clock 258, memory 260, 

digital terrain elevation data memory 262, and network interface 266. The clock may supply a 

plurality of timing signals for use by the server transceiver 252, the controller 256, the memory 

260, and digital terrain elevation data memory 262. The controller 256 is connected to the 

memory 260 and digital terrain elevation data memory 262 by a data bus 264. The controller is 

also coupled to a network interface 266 that enables the server 250 to communication with a 

larger network. The network may be a PCS network, cellular network, PSTN network, 

Bluetooth network, or other known wired or wireless type of network. 

[019] The instructions that are executed by the controller 256 and control the operation of 
the server 250 may be contained in memory 260. The memory 260 may be any type of memory 
or storage device that accepts and stores machine-readable instructions. The memory 260 and 
digital terrain elevation data memory 262 are shown as being separate memory, but in other 
implementations, only one memory may be employed. 

[020] A request for processing with the use of digital terrain elevation data is received from 
the SATPS receiver 102 by the server transceiver 252 via antenna 254 and processed by the 
controller 256 under the control of machine-readable instructions in memory 260. The request 
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from the SATPS receiver contains the code phases that are associated with the appropriate 
SATPS satellites 106, 108 and 110. The sever 250 may use the code phases from the SATPS 
receiver 102 to determine an estimated location. In other implementations, the request from the 
SATPS receiver 102 may contain an estimate of the location of the SATPS receiver 102. The 
controller 256 accesses the digital terrain elevation data memory 262 via data bus 264 to retrieve 
the digital terrain elevation data associated with the SATPS receiver's 102 estimated location. 
Examples of estimating the location of the SATPS receiver 102 at the server 250 include using 
either last known position of the SATPS receiver 102 or a triangulation scheme with directional 
antennas and known power setting of the radio transceiver 202. 

[021] The digital terrain elevation data from the digital terrain elevation data memory 262 
is then processed by controller 256 to determine the location of the receiver 102. Upon 
determination of the SATPS receiver 102, a message containing that location is sent to the 
SATPS receiver 102 by the server transceiver 252 via antenna 254. In another implementation, 
the digital terrain elevation data from the digital terrain elevation data memory 262 is sent to the 
server transceiver 252 by controller 256 and via antenna 254 and then transmitted to the SATPS 
receiver 102 and the navigation processor 208 determines the location of the SATPS receiver 
102 using the received digital terrain elevation data. In yet other implementations, the digital 
terrain elevation data may be contained at a common place within the network and the server 
250 would access the common place via network interface 266. 

[022] Examples of digital terrain elevation data include the 1999 NIMA standard digital 
dataset (DTED) level 0 for commercial and public use. THE DTED provides a worldwide 
coverage and is a uniform matrix of terrain elevation values that provide basic quantitative data 
of terrain elevation, slope, and/or surface roughness information. The DTED level 0 elevation 
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post spacing is 30 arc seconds (approximately 1 kilometer). In addition to the discrete elevation 

file, separate binary files may provide the minimum, maximum, and mean elevation values 

computed in 30 arc second square areas. Another example of digital terrain elevation data is 

GTOPO30 with 30 arc seconds spacing. The digital elevation model was derived from several 

raster and vector sources of topographic information. 

[023] The memory 212, 260 and 262 may be RAM, DRAM, SDRAM, EEPROM, a 
combination of RAM, DRAM, SDRAM, and EEPROM, or any other type of read/writeable 
memory. One skilled in the art will appreciate that all or part of systems and methods consistent 
with the present invention may be stored on or read from other machine-readable media, for 
example, secondary storage devices such as hard disks, floppy disks, and CD-ROMs; a signal 
received from a network; or other forms of ROM or RAM either currently known or later 
developed. Further, although specific components of SATPS 100 are described, one skilled in 
the art will appreciate that a positioning system suitable for use with methods, systems, and 
articles of manufacture consistent with the present invention may contain additional or different 
components. For example, the navigation process 208 may be a microprocessor, 
microcontroller, application specific integrated circuit ("ASIC"), discrete or a combination of 
other types of circuits acting as a central processing unit. 

[024] The digital terrain elevation data is processed along with the code phases to 
determine the location of the SATPS receiver 102. The navigation processor 208 forms four 
simultaneous equations with data retrieved from three SATPS spread spectrum signals and the 
digital terrain elevation data. The fourth equation derived from the digital terrain elevation data 
is a result of polynomial (in 2 variables of latitude <|) and longitude X) surface fit to the 
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appropriate terrain. To select the appropriate terrain from the digital terrain elevation data in the 
digital terrain elevation data memory 262 located in server 250, the SATPS spread spectrum 
signal 116, 118 and 120 from the three satellite 106, 108 and 110 are solved first for a fixed 
height "h" The fixed value of "h" may be initially assigned to the average value of "h" in the 
neighborhood of the base station (unlike known approaches of using the previous values of "h"). 
. Error in the fixed "h" is the absolute value of the difference between the average and the 
minimum or maximum value of "h". With this information, the three SATPS satellite position 
solution with fixed "h" comes with an estimated error ellipse. 

[025] Grid points along the direction of the major and minor axes of the error ellipse are 
constructed by the navigation processor 208. The step sizes (along the directions of major and 
minor axes) are made proportional (1.5 times) to the magnitudes of the major and minor axes 
respectively. In other implementations, other step sizes may be employed along the major and 
minor axes. The center of the ellipse is the two-dimensional location determined from the three 
SATPS satellite signals 116, 118, 120 and the fixed "h". In other implementations, the step 
size may be different. In the current implementation, 25 points are chosen in the grid. In other 
implementations, different number of points may be selected. Altitude values above the mean 
sea level at these points are obtained from the digital terrain elevation data (DTED) by indexing 
the four corner points in which the grid point resides and then using bilinear interpolation 
between the corner points. The obtained height values "H" are converted to World Geodetic 
Datum (WGS-84) "h" by adding the Geoid "N" separation derived from Earth Gravity Model 
(EGM-84) at the three SATPS satellite solution point determined from the three SATPS satellite 
signals 1 16, 1 18 and 120 and fixed "h". Earth Gravity Model EGM-96 or other types of earth 
gravity models may be employed for determining Geoid "N" separations. 
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[026] For the gird of 25 points of <|>, X and "h" a fourth order polynomial in <p and X is found 
using the Least Square (LSQ) method, resulting in 15 coefficients needing to be determined. To 
handle ill conditioning the polynomial is found in new variables that represent a scaled deviation 
from the center point (the three SATPS satellite solution determined from the three SATPS 
satellite signals 116, 118 and 120 and fixed "h"). A numerical method of Q-R decomposition is 
used, such as the modified Gram-Schmidt procedure, which makes Q only orthogonal rather 
than orthonormal (to avoid square root operations). 

[027] The maximum deviation of the grid point altitude from the surface fit is the error 
associated with this forth polynomial equation. If this error exceeds some error threshold, then 
the polynomial fit is rejected and more than one polynomial surface fit is required. In other 
implementations, different error thresholds may be employed. If the error threshold has not 
been exceeded, then only one polynomial surface fit is required. 

[028] The three polynomial equations associated with the three SATPS satellite signals 
116, 118 and 120, and the fourth polynomial equation obtained from the terrain elevation data 
are solved for curvilinear coordinates of <|>, X, H, and clock bias rather than using Earth Center 
Earth Fixed (ECEF) coordinates. The ECEF coordinate formulation is retained and changed 
from ECEF to WGS 84 curvilinear coordinates by the navigation processor 208 working with 
the corresponding Jacobian which is the coordinate transformation from differential ECEF 
coordinates to differential curvilinear coordinates. If there is a convergence then the converged 
solution is accepted if it belongs to the region of polynomial surface-fit. 

[029] Upon the location of the SATPS receiver 102 being determined and received by the 
SATPS receiver 102, the navigation processor 208, via interface 216 displays the location on 
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display 214. In other implementations, the digital terrain elevation memory 262 may be 
contained within digital terrain elevation data memory located within the SATPS receiver 102 as 
separate memory or may be within memory 212. In yet other implementations, the SATPS 
receiver 102 is only a module of a larger device, such as a cellular phone or other wireless 
device where the location data is output via location data output 218 to another module within 
the device. 

[030] In FIG. 3, a diagram of bilinear interpolation employed with digital terrain elevation 
data of FIG. 2 is shown. A table 300 of digital terrain elevation data of lxl degree increments is 
accessed from the South West Corner 302 (<f> r , X T ) by row 304 and column 306. Given a latitude, 
<|>u and longitude, X U9 the nearest South-West corner 302 of the table 300( DTED data file) is 
found and used as a reference to find an index in table 300. This index is used to retrieve the H. 
In other implementations, a different position or corner of the available digital terrain elevation 
data and a data type other than DTED may be used with appropriate changes to the equations. 
The equations are: 



[031] 



y 



u 




[032] 



where: 



y 



Reference Longitude for South West Corner of 
an available DTED data file. 

Apriori/User Longitude. 

Difference in degrees (longitude). 



[033] 




[034] 



where: 
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Reference Latitude for South West Corner of 
an available DTED data file. 

User Latitude. 

Difference in degrees (latitude) 



[035] 

[036] Where: 



brow = 



pacing 



y 

brow 



y * 3600 
A 2 

xl/ spacing 



DTED LEVEL 0 Grid Spacing of 30" Arc 
Seconds. 

Difference in degrees. 

Integer Row value within DTED data file in 
the range [0, 119]. 



[037] 

[038] Where: 



bcol - 



A<t>: 



'spacing 
X 

bcol 



3600 



A (J) 



spacing 



DTED LEVEL 0 Grid Spacing of 30" Arc 
Seconds. 

Difference in degrees. 

Integer Column value within DTED data file 
in the range [0, 119]. 



[039] The values of brow 304 and bcol 306 are used to find the index in the data file and 
then this index enables access to the altitude value H 3 , latitude and longitude X 3 . 
[040] Turning to FIG. 4 5 a diagram 400 of the four Orthometric heights identified in FIG. 3 
is shown. The Orthmetric heights Hi 408, H 2 406 and H 4 404 are obtained from (brow+\) and 
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bcol, (brow+\) and (bcol+\), and brow and (bcol+\) respectively. Then, H at <j> u and X u is 
obtained by: 



[041] 
[042] 



[043] 
[044] 



where: 



where: 



1=a Q> u - ^ 3 )3600 



Weighted Ratio from reference Latitude, Latitude 
grid spacing and DTED column location. 



y = 



(X u - i 3 )360Q 

spacing ~ br0W ) 



Weighted Ratio from user/reference Longitude, grid 
spacing and DTED column location. 



[045] 
[046] 



H^H X +{H 2 -H X )^ x'+(H 4 " #i ) • i + H 2 ~ H 2 - H 4 ) • x'*y' 

where: 

Hi,...,H 4 , represent 4 Orthometric heights in a given 
searched row and column result. 



Hi 



The interpolated Orthometric height and 25 of 
these points are determined. 



[047] To estimate the 25 points of "h" requires the estimation of the Geoid N from the 
EGM-84 model as a function of <|> u and X u , the processing is the same as determining H from the 
DTED. Once N is estimated, a linear shift by N is used to determine "h". 
[048] The LSQ solution from the three satellite signal equations (pseudoranges) and 
average "h" is determined by solving: 
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[049] P\ = p ]x -P x ) 2 + (s ly -p y ) 2 + (s lz - p z ) 2 • (1 - w, ) + b 
Pi =^(s 2x -P x ) 2 + (s 2y -p y ) 2 + (s 2z -p z ) 2 .(\-m 2 ) + b 

P3 = lj(hx -P X Y+ (hy -Pyf + (S 3s ~ P : f • 0 " ™ 3 ) + b 

h = ^P\-P X ) 2 +(p'y-Py) 2 +(P\~P;) 2 • Sgtl(h) 

[050] where (s ]x , s )y , s Jz ) are the ECEF coordinates of antenna phase center of a SATPS 
satellite "i" at the receive time, (p x , p y? p z ) are the ECEF coordinates of the antenna 1 14 phase 
center of the SATPS receiver 102, "b" is common offset in pseudorange measurements, m x is 
satellite mean motion correction term (given below), (p' x , p' y , p' z ) are ECEF coordinates of 
projection of (p x , p y , p z ) on the WGS-84 ellipsoid, p . is the measured pseudorange for SATPS 
satellite i, h is the height above WGS-84 ellipsoid and sgn(h)-l if h>0, sgn(h)=-l if h<0, and it 
is undefined when h=0 (here the equation itself reduces to an identity but the differential version 
of the equation is still defined). The "h" is given by the average "h" as obtained previously. 
The receive time is assumed to have errors of less than approximately 10ms so that the SATPS 
satellites 106, 108 and 1 10 positions as computed from the ephemeris have good accuracy. The 
mean motion correction term, m x is given as: 

[051] m i = - • (v, + co x s t ) o (s, - p) 

c 

[052] where x denotes vector cross product, ° denotes vector dot product, Vj is the velocity 
vector of SATPS satellite i, Sj is the position vector of SATPS satellite i, go is the Earth rotation 
vector and p is SATPS receiver 102 position vector, the x, y, and z coordinates of all vectors are 
in ECEF and all except a> correspond to the antenna 1 14 phase centers. 
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[053] If (p* x , p* y , p* z ) are the ECEF coordinates of the reference or approximate position 
which serves as the initial guess for (p x , p y , p 2 ) and let b* be the initial guess for pseudorange 
offset. Expanding the left hand side to a first order approximation using Taylor series gives the 
following equations: 



[054] 



[055] Where, lj is line of sight unit vector pointing from receiver to SATPS satellite i and d 
is down direction unit vector pointing along the downward normal to the WGS-84 ellipsoid, Ap x , 
Ap y and Ap z are differential position coordinates, Ab is differential pseudorange offset, Ap ] , 
Ap 2 and A/? 3 are differential pseudoranges, and Ah is differential h. The line of sight unit 
vector is given by: 











"Ap," 




"Ap," 


-h, 




-h, i 


• 


4p, 




Ap 2 


-h, 




-I* i 








&p 3 


-d x 


-<*, 


-d z 0 








Ah 



[056] 
[057] 

[058] 



the down direction unit vector is given by 

-cos A* »cos^* 
-sin X •cos^* 
-sin^* 



r* ( s , ~ P ) 



[059] where (<|>*, X*, h*) are the WGS-84 geodetic coordinates of p*. 
[060] Ap, = p, - p] 

[061] where p *\ are obtained from the right hand side of the equations of pseudoranges and 
average "h" and 
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[062] Ah = h - h* 

[063] The Taylor series equation is solved for Ap x , Ap y , Ap Zj and Ab. Then the estimates of 
position and clock bias are updated as: 



[064] 



Px 




~ * ™ 

P, 




"APx" 


Py 




p\ 


+ 


^Py 


Pz 




P. 




_ b _ 








Ab 



[065] This is the Newton-Raphson update. Next, the initial guess is replaced by the new 
estimate as: 



[066] 



* 

p, 




Px 


* 

p r 




Py 


* 

p, 




P.- 


b' 




b 



[067] and the iterations are continued until Ap x , Ap y , Ap z , and Ab become less than the 
threshold. 

[068] With the grid of 25 points a two-dimensional polynomial is set up in the auxiliary 
variables <f> and X which are given in terms of WGS-84 and X as: 

[069] ? = l=q<*-K) 

[070] where q is a scale factor (chosen as 100), and (|) c and Xc respectively are the WGS-84 
latitude and longitude in degrees of the solution obtained previously. The LSQ polynomial 
surface fit equation is given by: 

[071] h = p(<f>,Z) = c 0 -<f> +c, -X + c 2 -<j> 2 +<v^-X + c 4 -P +-~ + c m _ 2 - A" + c M _, 
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[072] The total number of coefficients, m, for degree, n are given by the recursive formula: 
m(n)=m(n-l)+(n+l) with m(0)=l. For degree, n=4, m=15. The coefficients, Cj, where i=1, 
m are obtained by solving the following linear equation using least square method. 



[073] 



<j> x A x tp 2 ^, • X x X 2 

02 A 2 (ft* <f> 2 ■ X\ 











c \ 




1 








K 




c 2 




1 








K 




• 












C4 






1 








A. 










. C «-l _ 





[074] where the subscript i (except on the coefficient) represents i-th point of the terrain. 
The points are chosen as follows: The center point (as given by Xc 9 and h c ) is the point that 
corresponds to the solution previously obtained. This solution also comes with horizontal error 
ellipse parameters of the semi-major axis, ae, the semi-minor axis, b e? and the angle the semi- 
major axis subtends with the east direction 0 C , measured anti-clockwise positive. This 
information is used to create a grid of points as: 





An 




COS0 e 


~s'mO e 




i • Aa e 


[075] 


Ae 




sm0 e 


COS0 e 


• 





i — Jc 5 * * * k j — k , * * ■ } k 



[076] 
[077] 

[078] 
[079] 



and 







<t> c 


• 











An/((N C + h)cos0 c ) 
Ae/(M c +h c ) 



where, 

M, 



a(\-e 2 ) 



(l-e 2 sin 2 ^) 



3/2 



and jV„ = 



*]\-e 2 sm 2 0 c 
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[080] where a is the semi-major axis of the WGS-84 ellipsoid and e is its eccentricity. Aae 
is chosen to be 1 .5 a e and Ab e is chosen to be 1 .5 b e • 

[081] The value of k is chosen to be 2, giving total number of points, r=25 and giving a 
coverage 3-sigma area along the semi-major and semi-minor axes. The degree of the 
polynomial n, is chosen as 4 th order giving number of coefficients, m=15. The system of 
equations therefore is 25 equations and 15 unknowns and is solved with the help of modified 
Gram-Schmidt procedure. 

[082] The LSQ polynomial surface fit equation in the usual matrix notation is A«C=H and 
the objective in least squares solution is to minimize (A«C-H) T »W-(A«C-H), where W is positive 
definite weighting matrix. The optimum solution is obtained by solving the set 
A T -W«A»C=A T# W«H. This set can be written as B T «B»C=B T »r«H, by using the decomposition 
W=r T «r and using B=A-r. This new set can further be written as R*C=D _, 'Q T -H, where B is 
decomposed as B=Q-R, with R unit upper triangular (diagonal elements of R are all ones and 
lower diagonals are all zeros) and such that Q T *Q=D, D being a diagonal matrix. The upper 
triangular set of equations can be solved using back-substitution method. In the above, two 
decompositions are used. The first is: W=r T «T. This can be done using Cholesky's method. 
Usually, W is diagonal and then so is T and it can be obtained simply by taking square roots of 
the diagonal elements of W. Even simpler case is when W=I, where I is identity matrix and then 
T=I as well. This equal weighting is used in the solution. The second decomposition is B=Q*R. 
This decomposition can be obtained by modified Gram-Schmidt method that gives Q, R and D 
by avoiding square root operations since Q is only orthogonal (not orthonormal). 
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[083] The equations from the three pseudoranges and polynomial surface equation are 

solved by LSQ solution method as described previously (using fixed "h") with two exceptions. 

The last (fourth) equation (fixed "h") is replaced by altitude equation as a polynomial in <|> and X, 

With this change, it is convenient to consider the first three equations in the unknowns of <|>, X, 

and h, rather than in the ECEF frame. So, the equations are written as: 



[084] />,=, 




-p y {<t>,A,h)) 2 


+ (% 


-p z W,A,h)) 2 
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-m t ) + b 






-p y (0,A,h)) 2 


+ ( S 2: 


-pM*M 2 


•0 


-m 2 ) + b 




l(s»-pM^h)) 2 + (s 3y 


-p y (t,A,h)) 2 


+ ( S 3: 


-pM^M 2 


•0- 


-m 3 ) + b 



/z-c 0 ^(^) + c i a(A) + c 2 ^(^) 2 +c 3 ^(^)aA + C 4 aW 2 +...-fc ffl _ 2 .A(X) n +c m _ x 
[085] where: 

Px = (N+h)»cos(()»cos^ 
p y = : (N+h)*cos())»sin)i 
p z =((l-e 2 )*N+h)»sin(() 

[086] The set of equations is solved using Newton-Raphson method. This time the initial 
guess is given by (c|) c , h c ) which is obtained previously. Further let b c be the initial guess for 
the pseudorange offset, taken again from the previous solution. Expanding the left side of the 
polynomial surface equations around the initial guess to a first order approximation using Taylor 
series and solving for A<|>, AX,, Ah, and Ab and then applying the same procedure as for the fixed 
"h" LSQ solution results in a sought solution. 

[087] In FIG. 5, a flow chart 500 of the steps for determining location with a SATPS 
receiver with digital terrain elevation data is shown. The SATPS receiver 102 starts 502 by 
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receiving three SATPS satellite spread spectrum signal 116, 118 and 120 and digital terrain 

elevation data from digital terrain elevation data memory 262 in server 250. Upon a reference 

location being determined, the orthometric heights at points 1 km apart in the easting and 

northing direction are retrieved with the reference location at the center 504. A total of (2-N+l) 2 

points are considered on a grid of (2-N+l)x(2*N+l) points. The orthometric "H" coordinates 

are converted to WGS 84 "h". The average "h" is determined and the "h" error equal to 

absolute value of the difference between the average "h" and the maximum or minimum "h" 

(whichever is the greater). 

[088] The SATPS pseudorange measurements, average "h" and the "h" error in step 506 
are used to find the position and corresponding horizontal error ellipse parameters (508). Then, 
with the position determined in step 508 at the center, "H" at points on a rectangular grid are 
constructed along the major and minor axes of the ellipse in 510. The interval between 
successive points along the major axis is 1.5 times the major-axis magnitude and along the 
minor axis, it is equal to 1 .5 times the minor-axis magnitude. A total of 25 points are considered 
on a grid of 5x5 points. 

[089] A two-dimensional polynomial of degree 4 is fit in the variables of § and X with 15 
coefficients to the 25 points previously obtained in step 510. The maximum residual error for 
the polynomial fit is determined and if the error exceeds a predetermined threshold in step 514, 
then an error message 516 is generated and processing is complete 518. Otherwise, if the 
predetermined threshold is not exceeded in step 514, then the SATPS pseudorange equations 
and the equation of the polynomial along with the maximum residual error of step 512 to 
determine the position 520. In step 522, a check is made to determine if the <|> and X obtained in 
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step 520 belong to the grid over which the surface-fit was done in step 512. If yes then the 

solution is accepted and processing is complete 518, otherwise the solution is rejected 524. 

[090] The processing has been shown as being stopped 518 for illustrative purposes. In 

practice, processing may continue and start at step 502 again, enter an idle state until instructed 

to process new location information, or execute other plurality of machine-readable instructions 

that are stored in some type of memory. The steps previously described, my occur with 

additional steps or with steps combined. 

[091] The foregoing description of an implementation has been presented for purposes of 
illustration and description. It is not exhaustive and does not limit the claimed inventions to the 
precise form disclosed. Modifications and variations are possible in light of the above 
description or may be acquired from practicing the invention. For example, the described 
implementation includes software but the invention may be implemented as a combination of 
hardware and software or in hardware alone. Note also that the implementation may vary 
between systems. The claims and their equivalents define the scope of the invention. 
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